library(plyr)
library(magrittr)
library(tidyverse)
library(readxl)


rm(list=ls())
home = 'C:/Users/Jason/Dropbox/VNA_Responsiveness/Analysis/JOP-dataverse/'


permutation.results = paste0(home, 'RI-analyses.Rds') %>%
  readRDS

perm.cit = subset(permutation.results, 
                  Result=='tab2col6' & term=='Citizen')$estimate
perm.firm = subset(permutation.results, 
                   Result=='tab2col6' & term=='Firm')$estimate
perm.cf.diff = perm.cit - perm.firm

perm.cit.speak = subset(permutation.results, 
                        Result=='tab3col7' & term=='Citizen')$estimate
perm.firm.speak = subset(permutation.results, 
                         Result=='tab3col7' & term=='Firm')$estimate
perm.cf.diff.speak = perm.cit.speak - perm.firm.speak



experimental.results = paste0(home, 'experimental-analyses.Rds') %>%
  readRDS

exp.cit = subset(experimental.results, 
                 Result=='tab2col6' & term=='Citizen')$estimate
exp.firm = subset(experimental.results, 
                  Result=='tab2col6' & term=='Firm')$estimate
exp.cf.diff = exp.cit - exp.firm

exp.cit.speak = subset(experimental.results, 
                       Result=='tab3col7' & term=='Citizen')$estimate
exp.firm.speak = subset(experimental.results, 
                        Result=='tab3col7' & term=='Firm')$estimate
exp.cf.diff.speak = exp.cit.speak - exp.firm.speak



dv_survey = paste0(home, 'survey_outcomes.xlsx') %>%
  read_xlsx %>%
  mutate(Treatment=factor(x=Treatment, 
                          levels=c('Control',
                                   'Citizen',
                                   'Firm')),
            Missing=as.integer(is.na(Q1))) %>%
  subset(!is.na(Treatment), 
         select=c(ID, 
                  Q1, 
                  Cultivation, 
                  Livestock, 
                  Treatment, 
                  FullTime, 
                  CentNom, 
                  Competitive))


educ = subset(dv_survey, 
              select=-c(Cultivation, 
                        Livestock)) %>%
  mutate(Education=1)
colnames(educ)[match('Q1', colnames(educ))] = 'Prepared'
cult = subset(dv_survey, 
              select=-c(Q1, 
                        Livestock)) %>%
  mutate(Education=0)
colnames(cult)[match('Cultivation', colnames(cult))] = 'Prepared'
live = subset(dv_survey, 
              select=-c(Cultivation, 
                        Q1)) %>%
  mutate(Education=0)
colnames(live)[match('Livestock', colnames(live))] = 'Prepared'
stacked = rbind(educ, 
                cult, 
                live) %>%
  subset(CentNom==0, 
         select=-CentNom)
colnames(stacked)[match('Treatment', colnames(stacked))] = 'Actual.Treatment'


assignments = paste0(home, 'RI_assignments.Rds') %>%
  readRDS %>%
  merge(stacked, 
        by='ID', 
        all=T) %>% 
  subset(!is.na(Actual.Treatment))


perm.cf.diff.edu = NULL
for(i in 1:1e4) {
  if(i %% 1e3 == 0) {
    cat('M\n')
  } else if(i %% 5e2 == 0) {
    cat('D\n')
  } else if(i %% 1e2 == 0) {
    cat('C\n')
  } else if(i %% 50 == 0) {
    cat('L')
  } else if(i %% 10 == 0) {
    cat('X')
  } else if(i %% 5 == 0) {
    cat('V')
  } else {
    cat('i')
  }
  tmpData = subset(assignments, 
                   Iteration==i)
  tmpLM = lm(Prepared ~ Citizen + Firm + Education*Citizen + Education*Firm + FullTime + Competitive, 
             data=tmpData)
  tmpDiff = tmpLM$coefficients['Citizen:Education'] - tmpLM$coefficients['Firm:Education']
  perm.cf.diff.edu = c(perm.cf.diff.edu, tmpDiff)
}; rm(list=c('tmpData','tmpLM','tmpDiff','i'))



stacked$Citizen = stacked$Firm = 0
stacked$Citizen[stacked$Actual.Treatment=='Citizen'] = 1
stacked$Firm[stacked$Actual.Treatment=='Firm'] = 1
actualLM = lm(Prepared ~ Citizen + Firm + Education*Citizen + Education*Firm + FullTime + Competitive, 
              data=stacked)
exp.cf.diff.edu = actualLM$coefficients['Citizen:Education'] - actualLM$coefficients['Firm:Education']



CI.cf.diff = CI.cf.diff.edu = CI.cf.diff.speak = NULL
for(i in seq(from=-2, to=2, by=0.001)) {
  CI.cf.diff = rbind(CI.cf.diff, 
                     c(shift=i,
                       shifted=exp.cf.diff+i,
                       upper=sum((exp.cf.diff+i)>=perm.cf.diff)/length(perm.cf.diff),
                       lower=sum((exp.cf.diff+i)<=perm.cf.diff)/length(perm.cf.diff)))
  CI.cf.diff.edu = rbind(CI.cf.diff.edu, 
                         c(shift=i,
                           shifted=exp.cf.diff.edu+i,
                           upper=sum((exp.cf.diff.edu+i)>=perm.cf.diff.edu)/length(perm.cf.diff.edu),
                           lower=sum((exp.cf.diff.edu+i)<=perm.cf.diff.edu)/length(perm.cf.diff.edu)))
  CI.cf.diff.speak = rbind(CI.cf.diff.speak, 
                           c(shift=i,
                             shifted=exp.cf.diff.speak+i,
                             upper=sum((exp.cf.diff.speak+i)>=perm.cf.diff.speak)/length(perm.cf.diff.speak),
                             lower=sum((exp.cf.diff.speak+i)<=perm.cf.diff.speak)/length(perm.cf.diff.speak)))
}; rm(i)
CI.cf.diff = as.data.frame(CI.cf.diff)
CI.cf.diff.edu = as.data.frame(CI.cf.diff.edu)
CI.cf.diff.speak = as.data.frame(CI.cf.diff.speak)



CIs.Primary = c((-1)*CI.cf.diff[CI.cf.diff$lower>0.025,][sum(CI.cf.diff$lower>0.025),]$shift, 
                exp.cf.diff,
                (-1)*CI.cf.diff[CI.cf.diff$upper>0.025,][1,]$shift)
CIs.Placebo = c((-1)*CI.cf.diff.edu[CI.cf.diff.edu$lower>0.025,][sum(CI.cf.diff.edu$lower>0.025),]$shift, 
                exp.cf.diff.edu,
                (-1)*CI.cf.diff.edu[CI.cf.diff.edu$upper>0.025,][1,]$shift)
CIs.Speak = c((-1)*CI.cf.diff.speak[CI.cf.diff.speak$lower>0.025,][sum(CI.cf.diff.speak$lower>0.025),]$shift, 
              exp.cf.diff.speak,
              (-1)*CI.cf.diff.speak[CI.cf.diff.speak$upper>0.025,][1,]$shift)



CIs = data.frame(rbind(CIs.Primary,
                       CIs.Placebo,
                       CIs.Speak)) %>%
  set_colnames(c('Lower',
                 'Estimate',
                 'Upper')) %>%
  mutate(Specification=factor(x=c('DV: Preparation\nNo central nominees + Covariates',
                                  'DV: Preparation\nNo central nominees + Covariates\nAccounting for Hawthorne effect',
                                  'DV: Speaking\nNo central nominees + Covariates'), 
                              levels=c('DV: Preparation\nNo central nominees + Covariates',
                                       'DV: Preparation\nNo central nominees + Covariates\nAccounting for Hawthorne effect',
                                       'DV: Speaking\nNo central nominees + Covariates')))
CIs

CIs.plot = ggplot(CIs, 
                  aes(x=Specification, y=Estimate, ymin=Lower, ymax=Upper)) + 
  geom_linerange(size=1) +
  geom_point(size=2) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  scale_y_continuous(breaks=seq(from=-0.3, 
                                to=0.3, 
                                by=0.1), 
                     labels=round(seq(from=-0.3, 
                                      to=0.3, 
                                      by=0.1), 
                                  digits=1)) +
  coord_cartesian(ylim=c(-0.3, 0.3)) +
  labs(x=NULL, y=expression(beta[Citizen]-beta[Firm])) +
  theme_minimal() +
  theme(panel.border=element_rect(color='black', fill=NA))


ggsave(filename='appendix-figure-14-1.png', plot=CIs.plot, path=home, width=9, height=3, units='in')
ggsave(filename='appendix-figure-14-1.tiff', plot=CIs.plot, path=home, width=9, height=3, units='in')
